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We develop a method to reconstruct the primordial power spectrum, P{k), using both temperature 
and polarisation data from the joint analysis of a number of Cosmic Microwave Background (CMB) 
observations. The method is an extension of the Richardson-Lucy algorithm, first applied in this 
context by Shafieloo & Souradeep l] . We show how the inclusion of polarisation measurements can 
decrease the uncertainty in the reconstructed power spectrum. In particular, the polarisation data 
can constrain oscillations in the spectrum more effectively than total intensity only measurements. 
We apply the estimator to a compilation of current CMB results. The reconstructed spectrum is 
consistent with the best-fit power spectrum although we find evidence for a 'dip' in the power on 
scales k « 0.002 Mpc~^. This feature appears to be associated with the WMAP power in the region 
18 < f < 26 which is consistently below best-fit models. We also forecast the reconstruction for a 
simulated, Planck-like Q survey including sample variance limited polarisation data. 
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I. INTRODUCTION 

With increasingly precise measurements being made of 
the CMB (Cosmic Microwave Background) or LSS (Large 
Scale Structure) it becomes progressively more important 
to determine how new observations could yield the great- 
est insight into processes occurring in the early universe. 
One such observable is the primordial power spectrum of 
curvature perturbations <i>(fc), 



t.3 ^ ^ 
P(fc) = —5^{k - fc')(*(fc)$*(fc')), 



(1) 



I a a i 

or other ex- 
\\m can mod- 



where k is the wavenumber. 

A generic prediction of the simplest inflationary mod- 
els is that the density perturbations should be adiabatic 
and near scale invariant. In these models the spectrum 
takes the form of a power law cx Current limits on 

this parametrisation place the spectral index rig ~ 0.963 
0. More complex inflationary models_ such as those 
with features on the potential [5, 
a small number of e-folds [ill. Tl2 
otic inflationary models [fl llSl. Tl6 
ify P(k) in a manner not compatible with the simple 
power law description. There have been many paramet- 
ric searches for the features pro duced by these models 
0, [H, 113, HH, m m, [13, m, M m m , although none 
have proved conclusive. On the other hand, there have 
been tantalising hints of anomalous features in the data, 
for example after the first year WMAP results were re- 
leased, there were strong indications of a cut-off in P{k) 
on large scales. With subsequent data releases the signif- 
icance of this feature has been reduced, although future 
observations of the polarisation of the CMB may provide 
more conclusive evidence Il3ll. 



A more thorough search for features in P{k) would be 
one in which the theoretical model is independent from 
the reconstruction. Methods such as the Richardson- 
Lucy deconvolution [ij, 123 . 1301 1. cosmic inversion 
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34 . 13511 ■ and other non-parametric approaches 
33, [HE m, EH, 113, Si all attempt to over- 
come theoretical bias. With current computing power 
these techniques are generally limited to recovering the 
spectrum for one set of assumed cosmological parameters. 
This allows one to use a fiducial CMB photon transfer 
function to integrate the primordial curvature perturba- 
tion into today's photon distribution perturbation. How- 
ever this assumption hides a significant degeneracy be- 
tween features in the primordial power spectrum and the 
physical parameters which determine the height and po- 
sition of acoustic peaks in the CMB when using only total 
intensity data in the reconstruction process [44] . As such 
it is not clear what the significance of any features found 
in the reconstructed P{k) should be. Adding polarisa- 
tion information into the reconstruction or inversion will 
significantly reduce this degeneracy since the response of 
polarisation CMB spectra is phase-shifted with respect 
to the total intensity response. 

The total intensity (T modes) of the CMB have been 
measured accurately by several instruments [sl. [isl [46ll47l. 
lisj to arcminute scales. Measurements in the gradient 
(i^-mode) and curl like (B-mode) polarisation compo- 
nents, and their correlation with the total intensity TE, 
and TB, lag behind due to their lower amplitude. How- 
ever a number of experiments are now measuring i?-mode 
polarisation with increasing signal-to-noise, starting with 
the first detection of EE spectrum [49l] and subsequent 
measurements H m m, [Hi [H. -B-mode polarisation 
has yet to be measured. Primordial tensor fluctuations 
may have an impact on the reconstruction of P{k) on 
large scales. However, _B-modes have not been detected 
yet and we neglect their contribution in this work. 

We show in this paper that when information con- 
tained in both total intensity and polarisation radiation 
transfer functions is used in the reconstruction of the 
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primordial power spectrum tighter constraints can be ob- 
tained. The paper is organised as foUows; in section HIl 
we introduce the extension of the Richardson-Lucy algo- 
rithm used to estimate P(fc). In Section Fill Al we test 
the algorithm using simple forecasts of CMB data. Our 
tests include template input models with radically bro- 
ken scale invariance. We explore current limits on P{k) 
in section IIIIBI and conclude in section |IVl 



II. 



AN EXTENDED RICHARDSON-LUCY Pk 
ESTIMATOR 



Direct primordial power spectrum reconstruction re- 
quires the inversion of the following relations 



XY 



dk 



Af(fc)Aj(fc)P(fc), 



(2) 



where X and Y represent T, E, or i3-type anisotropics, 
Cf^^ are the angular power spectra for the XY combina- 
tion and the Af{k) are the photon perturbation transfer 
functions. The transfer functions are obtained by inte- 
grating the full Einstein-Boltzmann system of differential 
equations. These describe the evolution of perturbations 
in the photon distribution functions in the presence of 
gravity and other sources of stress-energy. The functions 
determine all of the structure in the anisotropy spectra 
which arises after the initial conditions are set. Most no- 
tably the Cf^ contain distinct peaks due to the acoustic 
oscillation of the tightly coupled photon-baryon fluid in 
gravitational potential wells at the time of last scattering. 
The aim of any inversion method is to distinguish such 
features from any structure in the initial perturbation 
spectrum. 

For a finite sampling of the wavenumber space k Eq. ([2|) 
can be recast as an operator acting on the primordial 
spectrum Pk 



Ce — ^ FikPk, 



with operator 



= AlnfcAf.A 



(3) 



(4) 



where A In A: are the logarithmic k intervals for the chosen 
sampling. 

A solution for Pk cannot be obtained from a direct in- 
version of the as it is numerically singular. This is 
due to the high level of degeneracy in the transfer func- 
tions relating the power at any wavenumber k to angu- 
lar multipoles i although the system can be inverted by 
binning or smoothing appropriately to reduce the degen- 
eracy |52| . 

An alternative approach is the iterative inversion em- 
ploying the Richardson-Lucy (RL) method [H, H^l for 
image reconstruction. The RL method has been widely 



used in enhancing telescope images [55|, l56|, |57| and it can 
be shown that the RL method converges to the maximum 
likelihood estimator in the case of a Poisson distributed 
signal [5^ . In the following wc outline the use of the RL 
estimator in reconstructing the primordial power spec- 
trum as introduced by Shafieloo et al. Ij and extend it 
to include properly weighted contributions from polari- 
sation measurements. 

Consider the case where the original source plane is 
the isotropic, primordial Fourier space spanned by the 
wavenumber k and the convolved image plane is the space 
of angular multipoles £. In this case the convolution filter 
is Fik which relates modes the source power Pk to the 
image power Ci. 

The RL method provides an iterative solution to ^ 
for Pk, given an observed C^^^ with 



pi+l _ pi 



E 



/^obs 

F 



(5) 



where C| is the image obtained from the 1^^ iteration 
and Fik — Fik / J2e ^f-k such that in the hmit C} Cf"^ 
we have {Pl^^ — PI) ~> 0. It is also convenient to recast 
dni) as change in the Pk [l| 



^k - ^k 



^obs 

Fik ^ 



CI 



CI 



(6) 



such that the cut-offs in multipole space ^min and ^max 
are not propagated to the iterated solution through the 
broadening action of Fgk- 

As it is © does not account for errors in the observed 
image as it applies uniform weighting to all Cf°^. Since 
the RL estimator is not a well defined Maximum Likeli- 
hood estimator there is not a single choice of optimal 
weighting. Observational noise must therefore be in- 
cluded through an empirically chosen weighting function. 
We chose to add a weighting 



FikCl 



Gk = 



Y^Fik [Cl + at] 



(7) 



where cr^ is the reported error in the Cf°^. The weights 
Gk have the properties that if ai ^ C| then G/c — > 1 and 
if (T£ ^ Cj then Gk 0. The weighted estimator ([6]) 
now becomes 



pi+i 

^k 



PI 



/^obs /^i 

l + GkJ2Fik^' 



CI 



(8) 



Most CMB experiments observe only part of the sky 
which leads to correlated C^. The estimator does not 
account for the correlations as it includes only diagonal 
estimates of the uncertainty in Cf°^ . Indeed even full- 
sky experiments such as WMAP have correlated due 
to galaxy and source cuts. Experiments observing only 
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small fractions of the sky usually report results in as a set 
of bandpowers Cb = J2e ^e^e, where are the band- 
power window functions. The bandpowers are less cor- 
related than individual multipoles and allow us to apply 
the RL estimator to cut-sky experiments by calculating 
band filtered transfer function operators 



ek 



Fbk Pk 



(9) 



Using this we can define a generalised iterative estima- 
tor with contributions from any number of bandpowers b 
as 



^k 



PI 



1 + G 



'obs 



b 



bk- 



ci 



(10) 



where the W« 



5u' in the full-sky limit (WMAP) and 
are the reported bandpower window functions for other 
experiments. In this work we will be using bandpowers 
from ACBAR [H, QUaD BOOMERanG ^ and 
CBI [4^ which increase the range of scales probed past 
WMAP's resolution limits. 

A problem faced by the RL estimator when includ- 
ing polarisation measurements is that the polarisation 
bandpowers are correlated with the total intensity ones 
and cannot be included as further linear contributions to 
the sum in (|10p . Following the standard procedure for 
extending Maximum Likelihood estimators to polarised 
bandpowers we treat each multipole measurement as a 
2x2 matrix of bandpowers with 



(11) 



and similarly elevated the transfer function operator to 
a 2 X 2 matrix 



' bk 



E 



W^EEpEE 



(12) 



and we have disregarded any i?-mode contribution since 
current observational limits do not warrant its inclusion. 
We can now extend the RL estimator to include the full 
polarisation information as 



p«+i 
^k 



PI 1 



Tr 



b 



" bk 



(13) 



where is obtained by modifying ([7]) to sum over 
banded quantities and Fbk is normalised so that the sum 
over its diagonal elements is unity. 

We have implemented the estimator to reconstruct 
Pfc over a discretely sampled grid in k from a combination 



of CMB measurements listed above and for forecasted 
future data as described below. The full iteration typ- 
ically takes 0(10) minutes on a current desktop CPU. 
As a convergence criterion we iterate until the fractional 
change in the Pk solution is less than 0.01. Typically 
this takes 0(100) iterations to achieve. We reconstruct 
Pk over some 3000 k points between 6.9 x 10~^ Mpc~^ 
and 0.55 Mpc~^. The sampling is initially logarithmic 
and switches to linear at fc = 1.25 x 10~^ Mpc~^ follow- 
ing the typical setup in Einstein-Boltzmann solvers such 
as CAMB [53| but with higher resolution. 

As shown in [H the RL method requires a smooth- 
ing of the recovered Pk- We have found that applying 
a smoothing kernel at each iteration improves the con- 
vergence when using unbinned such as the WMAP 
results in the estimate. The smoothing helps to damp 
down large fluctuations at each iteration driven by the 
scatter in the observed data, due to either sample or noise 
variance. The behaviour is due to the estimator not be- 
ing a well defined maximum likelihood one and thus not 
taking into account the full correlation structure of the 
data in its weighting. We have chosen a simple smooth- 
ing kernel for the applications reported in this work by 
taking a nearest neighbour running average over the Pk 
between each iteration. Given our k sampling this de- 
fines a smoothing size of A log fc « 0.04 Mpc"i at low k 
and Afc « 2 x 10"** Mpc~^ at high fc. Other smooth- 
ing kernels may lead to more efficient behaviour of the 
estimator such as faster convergence but would be more 
computationally demanding and would lead to a further 
loss of resolution in k space. 

Since the estimator requires 0(10) minutes to converge 
it is not currently possible to employ it as part of a full 
Markov Chain Monte Carlo parameter search in place of 
the standard primordial spectrum parameter; amplitude 
As and spectral tilt Us ■ This may be possible in future if 
the algorithm is optimised and/or parallelised but in this 
work we chose to explore the Pk "source" plane assuming 
a set of fixed fiducial transfer functions based on the best- 
fit parameters. The transfer functions are obtained from 
CAMB with parameters flbh'^ = 0.0226, l^^/i^ = 0.108, 
e = 1.041, and T = 0.076. 

To estimate the confidence limits around the sampled 
Pk solution we Monte Carlo the iterative estimator by 
generating 1000 simulated sets of C^^^. For simplicity we 
assume Gaussian distributions for all the C^^'^ since more 
accurate distributions are not easily reconstructed from 
some of the published data. We iterate to convergence on 
each of the simulated data sets and obtain a covariance 
matrix for the samples by averaging over the ensemble 



n=N 



'kk' 



N 



1 



{Pk)W - (Pk')), (14) 



with {Pk) = 'l2n=i Pk/N. We quote 1 and 2-a errors 
by taking multiples of the square root of the diagonal 
elements of the cr^j,, matrix. 
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FIG. 1: The reconstruction of several test spectra. The test models used to generate the simulated Ce, shown in black (solid) 
curves, are (a) A 10% decrease in power from the WMAP5 best fit amplitude, (b) the WMAP5 best fit model including 
running dris/dhik = —0.037, (c) a localised feature at around k = 0.02 Mpc~^, and (d) a model with sinusoidal oscillations 
superimposed on the best fit power law spectrum. The green (solid) curves show the best fit spectrum used as initial guess 
in the iteration. The red (dotted) curves are the converged reconstructions using only total intensity data whereas the blue 
(dashed) curves use both total intensity and polarisation data. The Ce forecasts assumed an experiment with similar properties 
to Planck. 



III. APPLICATION OF THE EXTENDED RL 
ESTIMATOR 

A. Tests on forecasted CMB data 

We have tested the extended RL estimator on a num- 
ber of input Pk templates using forecasted data with 
similar experimental properties as the upcoming Planck 
satellite mission 0. We assume a total of 12 detectors 
with NETs of GifiK/^ observing 80% of the sky over 
12 months with a resolution of 7 arcminutes FWHM. We 
calculate errors around our fiducial CMB best-fit mod- 
els in both total intensity and polarisation spectra for 
this experimental setup and use these together with 
samples on the fiducial model to test the estimator's con- 
vergence. We consider multipoles of ^ < 1900 for total 
intensity spectra and £ < 1000 for polarisation. We have 
not taken into account any residual error from foreground 
subtraction in our forecasts. Thus our forecast are on the 



optimistic side of the accuracy achievable, particularly in 
the polarisation where foreground removal will certainly 
have a significant impact on errors at ^ < 1000. 

We use four separate Pk templates to generate the Cgs. 
The first is a standard power law with Us — 0.96 but 
with a lower amplitude than the best-fit, the second is 
a running model with dns/dlnk = —0.037, the third 
is a power law with a sharp, compensated feature at fc = 
0.02 Mpc-i 0, 0, B M, M, M,m, [13 and the fourth 
is apower law with superimposed sinusoidal oscillations 
0, M, IHj [s^j [Ml • The extended estimator is run on all 
four sets of Ci forecasts and the resulting Pk solution is 
compared to the input one. In Fig. [T]we show the results 
for runs including Cf^ only and runs including Cf^ , 
Cf^, and Cf ^. 

Each run is started with a first guess Pk (green/long- 
dashed) which is the current best-fit power law spectrum. 
In all cases the structure in the input Pk is reproduced to 
some degree over a large range in k. Beyond the k range 
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FIG. 2: Confidence regions around tiie reconstructed Pk for tlie same test models sliown in Fig. [T] Tfie sfiaded areas indicate 
tlie 1 and 2-a confidence regions obtained from 1000 Monte Carlo realisations of the observed Ci. The red (solid contours) 
shows the result for polarisation data included overlaid on the blue (dashed contours) showing the result for only TT included. 
The inclusion of polarisation has the largest impact in the case with superimposed oscillations. 



shown in Fig. [T] there are significant departures from the 
correct solution. On the large scales (fc < 0.0001 Mpc~^) 
this is due to the cut-off in the transfer functions due 
to the horizon scale. At small scales {k > 0.1 Mpc~^) 
the cut-off is due to the resolution limit of the forecasts. 
We find the method is particularly well suited for recon- 
structing long wavelength structure such as in the run- 
ning and oscillating model. Although the sharp feature in 
the third model is present in the reconstructed spectrum, 
its amplitude is not recovered accurately. This example 
highlights the limitations of such methods in reconstruct- 
ing high frequency features in the primordial spectrum. 
In the first model we also see high frequency features on 
small scales, this is due to the difference in amplitude be- 
tween the initial guess and correct solution. A power law 
initial guess with closer starting amplitude to the correct 
value reduces this small scale noise. This could be easily 
obtained by carrying out a standard power law fit to the 
data before running the reconstruction estimate. 

In Fig. [2] we show the same set of reconstructions but 
include 1 and 2-a confidence regions obtained from the 
Monte Carlo covariance matrix. We have over sampled 



the Pk and care should be taken in interpreting the signif- 
icance of any feature as the samples are highly correlated. 
The plot includes a set of contours for the reconstruc- 
tion using only total intensity data (blue/dashed) and 
including polarisation (red/solid). The range in k has 
been modified to emphasise the regions where the best 
limits are obtained. The addition of polarisation data 
gives additional constraints on the smaller scales partic- 
ularly in the oscillating model case. This is not surprising 
since total intensity and polarisation measurements are 
complementary in differentiating any structure in the pri- 
mordial spectrum from acoustic oscillations imprinted on 
the CMB at last scattering. The polarisation data has 
little effect in improving the constraints on the spectrum 
with a compensated feature although a detection of the 
feature is evident. 



B. Current Limits 

In this section we use the extended RL estimator to 
reconstruct the primordial power spectrum using cur- 
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FIG. 3: Current limits from a combination of CMB data sets 
(WMAP, ACBAR, QUaD, BOOMERanG and CBI). There is 
some evidence of a dip in power at around k « 0.002 below 
the best fit power law model. Shaded regions are defined as 
in Fig. [2] 



FIG. 4: As an indication of the origin of the dip at fe « 0.002 
Mpc^^ we remove the data between ^ = 18 and £ — 26 and 
re-run the estimator. The red/yellow contours show the effect 
of the removal over the original estimate (blue/cyan). 



rently available data. Our combination of data includes 
the latest WMAP release ^ including TT for ^ < 700, 
the CBIpol results |48i] including TT band powers 2 to 
12 and all polarisation band powers, the BOOMERanG 
2003 [iTi flight band powers, the QUaD ^ results, and 
the ACBAR 2008 [45] excluding the last two band powers 
due to the excess power. We have excluded the WMAP 
polarisation measurements as they are too noisy to be 
used in the reconstruction would require binning to re- 
duce the scatter. 

The overall effect of adding in the sub-orbital exper- 
iment is to extend the k range of the reconstruction 
past k 0.07 Mpc"i where the WMAP noise becomes 
too large for the data to be included in the estima- 
tor. Given WMAP's resolution in multipole space, it 
strongly dominates the reconstruction on scales larger 
than k ~ 0.07 Mpc~^. In Fig. [3] we show the results of 
the reconstruction. There is no indication of a departure 
from a pure power law for the recovered spectrum except 
for a dip at fc ~ 0.002 Mpc"^ 

To determine the origin of the feature at fc '--^ 
0.002 Mpc~^ we re-estimate the power spectrum with 
a cutout of the WMAP data in the range 18 < < 26 
where the Cg lie significantly below the best-fit models. 
Fig. |4] shows the result of for the reconstruction. No ev- 
idence for the feature remains when the 18 < £ < 26 
WMAP data is removed from the estimate indicating 
that the dip is associated with the feature in the Ce. 



IV. DISCUSSION 

The RL estimator provides a model independent 
method for reconstruction of the primordial power spec- 
trum of perturbations Pk from measurements of the CMB 
angular power spectrum C^. We have extended the RL 
estimator to use on multiple data sets including properly 
weighted polarisation data and have used Monte Carlo 
realisations of the input Ci measurements to estimate 
confidence limits around the reconstructed spectra. 

We have applied the new estimator to current measure- 
ments. These include band powers from the QUaD experi- 
ment which provide the highest signal-to-noise measure- 
ments of EE power so far. Including the polarisation 
information increases the constraints on the power spec- 
trum reconstruction as it carries independent information 
which is complementary in phase to that of the total in- 
tensity. However current TE and EE measurements are 
still noise dominated and do not contribute significantly 
but future sample variance limited measurements will 
help to constrain any structure in the primordial spec- 
trum as shown in our examples using forecasted data. 
An exception to this is the current WMAP polarisation 
data which may have some impact if it were binned to 
reduce noise scatter. We have not explored this option in 
this work as this would have required binning of the TT 
data too with a resulting loss of resolution in k space. 

We have found that the addition of polarisation data is 
particularly helpful in constraining oscillatory structure 
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FIG. 5: We show the improvement that polarisation data can 
have on the estimation of the spectrum. The blue area is the 
uncertainty on given only the WMAP 5-year data. In the 
foreground we have combined this data with polarisation data 
from our Planck-like experiment and show that it significantly 
reduces the inherent error. 



in Pk- Spectra with superimposed sinusoidal features 
have been considered in the hterature and have been con- 
strained using model dependent fits d, [6l|, [gI, [G^I . 

Due to the empirical nature of the weighting used the 
estimator has a limited acceptance range in the signal- 
to-noise of the band powers it is run on. Low signal-to- 
noise Cg^'^ or C^^^ cause the estimator to converge very 
slowly. Conversely if the data is weighted too strongly 
there is a loss of resolution and signal-to-noise in the re- 
constructed Pk ■ The best approach is to exclude the low 
signal-to-noise measurements from the data being used. 



It is possible that other weighting schemes may be more 
efficient in allowing the whole range of Ci measurements 
to be used. However, in its present form, the estimator 
will be most useful when all polarisations will be sample 
variance limited out to the same £. As an example of the 
improvement this will have over current estimates we add 
our forecasted TE and EE data to the current measure- 
ments and run the extended RL estimator. The result is 
shown in Fig. [5] which shows an overall improvement in 
the Monte Carlo confidence limits over the entire range 
in k being reconstructed. 

In this work we have not explored the degeneracy of 
the method with the physical parameters determining 
the structure of the CMB transfer functions. The cur- 
rent run-time to convergence for each reconstruction is 
too long to allow the tens of thousands of runs required 
to implement the method as part of a full Markov Chain 
Monte Carlo exploration of the entire parameter space. 
However with parallelisation and more efficient conver- 
gence it may become possible to do this in the near fu- 
ture. 

There are further extensions of the RL estimator that 
could increase its effectiveness. The addition of other ob- 
servables with different transfer functions such as galaxy 
redshift surveys or cosmic shear surveys will provide com- 
plementary information in the reconstruction. Our ex- 
tension of the RL estimator is compatible with these 
other observables as they can be included as further inde- 
pendent data points with properly formatted and binned 
transfer functions. Such extensions would contribute key 
information in the reconstruction on scales where the 
presence of foregrounds in the CMB reduces its effec- 
tiveness. 
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